﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace SmartMathLibrary
{
    /// <summary>
    /// This class represents a Moore-Penrose inverse matrix for real numbers.
    /// </summary>
    [Serializable]
    public class MoorePenroseInverseMatrix : Matrix
    {
        /// <summary>
        /// Initializes a new instance of the <see cref="MoorePenroseInverseMatrix"/> class.
        /// </summary>
        /// <param name="rows">The number of rows of the matrix.</param>
        /// <param name="columns">The number of columns of the matrix.</param>
        public MoorePenroseInverseMatrix(int rows, int columns)
            : base(rows, columns)
        {
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="MoorePenroseInverseMatrix"/> class.
        /// </summary>
        /// <param name="matrix">The matrix to clone.</param>
        public MoorePenroseInverseMatrix(Matrix matrix)
            : base(matrix)
        {
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="MoorePenroseInverseMatrix"/> class.
        /// </summary>
        /// <param name="matrix">The matrix to clone.</param>
        public MoorePenroseInverseMatrix(MoorePenroseInverseMatrix matrix)
            : base(matrix)
        {
        }

        /// <summary>
        /// Calculates the Moore-Penrose inverse matrix of the specified matrix.
        /// </summary>
        /// <param name="matrix">The matrix from which the Moore-Penrose inverse matrix should be compute.</param>
        /// <returns>The computed Moore-Penrose inverse matrix.</returns>
        public static MoorePenroseInverseMatrix FromMatrix(Matrix matrix)
        {
            if (matrix == (Matrix) null)
            {
                throw new ArgumentNullException("matrix");
            }

            SingularValueDecomposition decomposition = new SingularValueDecomposition(matrix);

            decomposition.ExecuteDecomposition();

            double[] diagonal = new double[decomposition.S.Count];

            for (int i = 0; i < decomposition.S.Count; i++)
            {
                if (decomposition.S[i] != 0)
                {
                    diagonal[i] = 1 / decomposition.S[i];
                }
            }

            Matrix diag = new Matrix(decomposition.S.Count, decomposition.S.Count);

            for (int i = 0; i < diagonal.Length; i++)
            {
                diag.SetValueAtPosition(i, i, diagonal[i]);
            }

            Matrix u;
            Matrix v;

            if (decomposition.Transposed)
            {
                u = decomposition.U;
                v = decomposition.V;
            }
            else
            {
                u = decomposition.V;
                v = decomposition.U;
            }

            Matrix result = u * (diag * v.Transpose());

            if (decomposition.Transposed)
            {
                result = result.Transpose();
            }

            return new MoorePenroseInverseMatrix(result);
        }
    }
}